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We introduce the sandpile model on multiplex networks with more than one type of edge and 
investigate its scaling and dynamical behaviors. We find that the introduction of multiplexity does 
not alter the scaling behavior of avalanche dynamics; the system is critical with an asymptotic power- 
law avalanche size distribution with an exponent r = 3/2 on duplex random networks. The detailed 
cascade dynamics, however, is affected by the multiplex coupling. For example, higher-degree nodes 
such as hubs in scale-free networks fail more often in the multiplex dynamics than in the simplex 
network counterpart in which different types of edges are simply aggregated. Our results suggest 
that multiplex modeling would be necessary in order to gain a better understanding of cascading 
failure phenomena of real- world multiplex complex systems, such as the global economic crisis. 
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I. INTRODUCTION 

Complex network theory has provided a useful tool 
for studying the structure and the dynamics of complex 
systems [TJ [2] . Various dynamical processes occuring on 
networks have also been extensively studied [3]. An im- 
portant class of dynamical processes is the so-called cas- 
cading failure dynamics, addressing how the failure of 
nodes induces other nodes' failure over the network. Var- 
ious real- world phenomena such as blackouts of electrical 
power grids [H |5], Internet congestion cascades [H |6], 
or global economic crises [7 can be modeled under the 
framework of cascading failure on networks. Prototypi- 
cal models of cascading failure include the sandpile model 
[8] , the threshold cascade model of behavioral adoptions 
[9], and the capacity-based overload cascade model [5 . 

Most studies thus far have focused on networks with 
a single type of link, the so-called simplex networks. In 
many real- world complex systems, nodes participate in 
more than one type of interaction with other nodes in 
the system; People in society have friendship, family ties, 
professional collaborative relations, etc. Countries in the 
global economy interact through various financial and po- 
litical channels; they are multiplex networks. Different 
interaction channels do not always operate in the same 
way, nor do they operate in a completely autonomous 
way. Therefore, the multiplexity of networks introduces 
another nontrivial structural factor in complex networks, 
the dynamic consequence of which is to be understood. 
In this perspective, there have been recent studies related 
to the multiplexity of networks in the form of interdepen- 
dency between network layers [10], interacting networks 
[11], and multi-relational structure of online social net- 
works [12]. These studies showed that the multiplexity 
presents nontrivial structural patterns and connectivity 
[TTHT3] and can drastically change the robustness prop- 
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erties of interdependent network systems [10] . This sug- 
gests that the simplex network framework would gener- 
ically be incomplete and thus one needs to consider the 
multiplex network framework for a better understanding 
of complex systems. 

Each type of interaction in multiplex networks defines 
a network layer; for example, a social network may com- 
prise a friendship layer, a family tie layer, etc. A change 
of a node's state may be determined by intralayer dy- 
namics; for example a person in a social network may 
adopt a behavior (using a smartphone app, say) by the 
peer pressure from his/her friends or by a demand due 
to the work environment alone. However, the potential 
influence of his/her state change is not confined within 
that layer; even though you adopted the app due to your 
work demand, your usage of the app may spread to your 
friends or family as you may recommend it through your 
friendship and family tie layers. In this way, the layers of 
a network get coupled and become interdependent. Such 
a multiplex coupling is generically not equal to a simple 
aggregation of link types, thereby leading to nontrivial 
dynamical consequences. 

In this paper, we apply multiplex modeling to a sim- 
ple model of cascading failure, the sandiple model, and 
investigate its dynamical impact. We begin by review- 
ing the sandpile dynamics on random, simplex networks 
in Sec. II. The multiplex sandpile model is introduced 
in Sec. Ill, and the principal results will be presented in 
Sec. IV. Other model variants are considered in Sec. V, 
and the paper concludes in Sec. VI. 



II. SANDPILES ON RANDOM NETWORKS 

The sandpile model was originally introduced by Bak, 
Tang, and Wiesenfeld (BTW) in 1987 [8]. Since then, this 
model has been an archetypical theoretical model for in- 
vestigating cascading failure dynamics and self-organized 
criticality [14 . The dynamics proceeds as follows, i) 
Each node is assigned a prescribed threshold height, usu- 
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ally taken to be its coordination number {cx = z^}. ii) 
At each step, a grain is added at a randomly chosen node, 
say X, raising the node's height by one, hx iii) 

If the height of the node becomes equal to or exceeds 
its threshold, that is hx > Cx, the node becomes un- 
stable and topples Cx grains to its neighbors, such that 
hx hx — Cx and hy ^ hy -\- 1 for y being the neighbor 
of X. iv) If this toppling event induces instability in some 
neighbors, they also topple in the same way; a cascade of 
toppling events (an avalanche) occurs until there remain 
no more unstable nodes in the system, v) Repeat ii)-iv). 
BTW showed that without a fine-tuning of the external 
parameters, this model could exhibit a scale invariance 
in the form of a power-law distribution of avalanche sizes 
(the number of nodes participating in a given avalanche) , 
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hence, the name self-organized criticality [8 . 

In Euclidean lattices in low dimensions, the BTW 
sandpile model exhibits nontrivial critical behaviors [14^ . 
Its dynamics on random networks, however, can be un- 
derstood more easily, following the mean- field behaviors. 
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[T6] . Later, the model has been con- 



sidered on the scale-free networks with a power-law de- 
gree distribution, P{k) ~ where the degree k is the 
number of links of a node pm2Q] . For the BTW sandpile 
model with q = ki on scale- free networks, the unconven- 
tional mean-field behavior was obtained such that the 
avalanche size exponent r is given by [17 
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(7 > 3), 
(2 < 7 < 3). 
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Here, we briefly present the analytical approach based 
on the branching process approximation leading to these 
results [171 IS, • Assuming that the avalanche proceeds 
without forming a loop, it can be mapped onto a branch- 
ing process. In this mapping, the probability to make 
a /c-branch {k > 1) at each branching event is given by 
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The first factor comes from the likelihood that a node 
with degree k would gain a grain from its unstable neigh- 
bor, which is proportional to its degree. The second fac- 
tor comes from the assumption that the timing at which 
the node gets a grain from its neighbor is random, such 
that the probability to become unstable (that is, h = k—1 
at the moment) is 1/k. qo is given by the normalization 
as go = 1 — YlV=i^ki and the mean branching number 
B = XlfeLo is obtained as 
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Therefore, the branching process is critical. The asymp- 
totic (large s) behavior of the avalanche size distribution 



P{s) can be obtained from the singular behavior near 
unity in the corresponding generating function [18 . On 
scale-free networks, the branching probability decays al- 
gebraically, qk ~ P{k) ~ k~^ ^ which gives rise to a singu- 
lar term ~ (1 — o;)^"^ in its generating function, leading 
to the asympototic power-law behavior of P{s) for 7 < 3 
as given by Eq. (2). Thus, the behavior of the branching 
probability qk provides an essential clue to the sandpile 
dynamics. 



III. SANDPILE MODEL ON MULTIPLEX 
NETWORKS 

In this paper we generalize the BTW model on mul- 
tiplex networks. In multiplex networks, each type of in- 
teraction functions as a network layer within which the 
dynamics takes place as if in the usual single network 
case. However, the state change of a node in one layer ap- 
plies also to other layers, giving rise to couplings between 
layers' dynamics. Such an interlayer coupling induces in- 
terdependency between layers and renders the multiplex 
dynamics rather nontrivial JOl EH [22] . To illustrate the 
effect of multiplexity in its simplest setting, we consider 
in the paper multiplex networks with two kinds of inter- 
actions or a duplex networks. 

Specifically, the BTW sandpile model on duplex net- 
works is run as follows: i) Each node x is assigned the 
thresholds c^^x equal to its degree for each layer; that is, 
= ki^x for = 1 and 2. ii) At each step, a grain is 
added to a randomly chosen node x in a randomly cho- 
sen layer ^; h^^x ^ h^^x + 1- ih) If the height of this 
node reaches or exceeds its threshold in any layer, that 
is, either hi x > ci x or h2 x ^ ^2x1 the node becomes 




FIG. 1: (Color) Schematic illustrations of the model dynam- 
ics, (a) A network is composed of two layers of red (top) and 
blue (bottom) interactions, respectively, (b) In multiplex dy- 
namics, the individuality of different interactions is respected; 
thus, each node has two piles of grains that follow the multi- 
plex 102 dynamics, (c) In the simplex 102 case, two types 
of interactions are aggregated, and the conventional network 
sandpile is run on the aggregated network. 
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unstable, initiating an avalanche. The layer in which the 
threshold is exceeded is called the unstable layer, iv) 
Each unstable node topples on both layers; in the unsta- 
ble layers, it topples ci^x grains, one for each neighbor, 
whereas it topples only h^'^x grains in the other layer, giv- 
ing one grain to each of a random subset of hi^^x neighbors 
[Fig. 1(b)]. That is, 

hi,x ^ hi^x - min[Q,^, h^^x] (1(^2 model), (5) 

and the height of all the nodes receiving a grain increases 
by one. v) If this toppling event causes other nodes to be 
unstable, these unstable nodes topple in the same way as 
in iv) . vi) Repeat iv)-v) until there are no more unstable 
nodes in the system, and the avalanche ends. During 
each toppling, the grain is lost with a small probabilty, 
taken to be / = 10~^ in this paper, which is necessary to 
prevent overloading the system. We call this multiplex 
dynamics model the 102 model. 

In our l(g)2 multiplex setting, a node's failure (top- 
pling) due to a specific layer dynamics induces its simul- 
taneous failure in other layers regardless of its conditions 
in those layers. This setting is motivated by a number 
of real-world situations such as the global economy, in 
which one country's recession due to an imbalance in the 
trade system can affect the country's failure in a credit 
network, for example. The 102 model is a stylized model 
of such scenarios, aiming to address the effect of multi- 
plexity. To highlight the effect of multiplexity in the 102 
model, we compare it to a non-multiplex, or simplex, 
counterpart of the model, which we call the 102 model. 
In the 102 model, we ignore the individuality of differ- 
ent interactions and simply overlay two layers [Fig. 1(c)] 
to form a simplex aggregate network, on which the con- 
ventional BTW model is run. In the next section, we 
study and compare the avalanche dynamics of these two 
models. 



IV. RESULTS 

The primary quantity of interest in sandpile dynamics 
is the avalanche size distribution P(s), giving the likeli- 
hood that an avalanche of size s is to occur in the steady 
state of the model [H [14] . We meaure the avalanche size 
in terms of the number of nodes becoming unstable dur- 
ing an avalanche. A related quantity is the number of 
toppling events during an avalanche, which scales sim- 
ilarly when the formation of loops during an avalanche 
is negligible, which is the case for sandpile dynamics on 
random graphs [17]. One can also measure the distribu- 
tion of avalanche durations, which also follows a power 
law, whose exponent is related to the avalanche size ex- 
ponent by a scaling relation In this work we focus 
exclusively on the avalanche size distribution. 
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FIG. 2: (Color online) Avalanche size distribution P{s) of 
the 102 (circle) and the 1©2 (square) model on duplex ER 
networks with N = 10^ and mean degree z = 4 per layer. 
The slope of the dotted line is —3/2 and is drawn for the eye. 
Data are logarithmically binned. 



A. On Duplex Erdos-Renyi Layers 

We first consider the sandpile dynamics on multi- 
plex networks with two Erdos-Renyi (ER) f23| layers of 
N = 10^ nodes and the equal mean degree z = 4. Each 
node has degrees ki and /c2 in the first and the second 
layers, respectively, each of which is a Poisson random 
variable with the mean equal to z = 4. To obtain the 
avalanche size distribution and related statistics, we sam- 
pled 10^ successive avalanches after the initial transient 
period. For the simplex 102 model, the substrate net- 
work is nothing but an ER network with mean degree 
(k) = 2z = S. Therefore, the 102 model is trivially crit- 
ical, with the standard mean- field exponent ri02 = 3/2 
[m [S] 5 which is also confirmed by numerical simulations 
(Fig. 2). 

For the multiplex 102 model, the sandpile dynamics 
would become more complicated due to layers' coupling 
or interdependency. Even the criticality of the model 
is not obvious a priori The numerical simulation of 
the 102 model, however, showed that its avalanche size 
distribution is essentially identical to that of the 102 
model; P{s) follows a power law, that is, the model is still 
critical, with the same mean- field exponent ri^2 =3/2 
(Fig. 2). Such a robust mean- field behavior was also re- 
ported recently for sandpile models on interdependent 
networks in somewhat different settings [24 . 

Although the two models show similar P{s) with the 
same exponent r, the specific sandpile dynamics un- 
derlying this distribution are different, given the pres- 
ence/absence of multiplexity. Compared to the simplex 
case that can be understood in terms of the framework 
presented in Sec. II, the multiplex dynamics follows a 
more complicated cascade process. To provide more de- 
tailed insight, we again employ the branching process ap- 
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FIG. 3: (Color online) Branching probability qk for the 102 
(circle) and the 102 (square) model. (Inset) Plots of the 
factor 0{ki) in Eq. (7) for the two models. The slopes of the 
dotted lines are —1.0 (upper) and —1.3 (lower), respectively 
and are drawn for the eye. 



proximation following Sec. II. For the simplex 102 case, 
the branching probability is given by, from Eq. (3), 



102 
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where Pi^2{k) is the degree distribution of the aggregate 
network of two layers. Equation. (6) is confirmed by nu- 
merical simulations [Fig. 3(a)]. For the multiplex 102 
model, the associated branching process can be divided 
into two parts: one through the unstable layer (denoted 
layer 1) and the other through layer 2. The latter is 
the consequence of multiplexity and the layers' interde- 
pendency. To obtain the branching probability into k 
branches, one has to consider the composite branching 
of the /Ca^-branch in layer 1 and the /i^^-branch in layer 2. 
Therefore, one can write 



ki = l 



ki<k 



k2 



The first summand gives the probability that a node with 
degree ki receives a grain and becomes unstable, making 
a /ci -branch. The second summation gives the probability 
that the node has k — ki grains in layer 2 at the moment 
(thereby k — ki branching), assuming again the random 
timing of topplings. Summing the combined factor over 
/ci, we have qk. Multiplexity is clearly manifested in the 
second factor, but it does play a role in the first factor, 
too. The unspecified factor 6{ki) addresses the probabil- 
ity that a node becomes unstable as it receives a grain, 
which was 1/k for simplex sandpile, Eq. (3). The as- 
sumption behind it was that all heights are equally likely 



FIG. 4: (Color) Plot of the normalized difference between 
the frequency of node failure as a function of node degrees, 
0(A:i,A;2), defined in the text. A positive (red) value of (j) 
means that the node fails more frequently in the 102 model 
than in the 102 model. A negative (blue) (f) indicates the 
opposite. 



when encountering a toppling, as the grains are succes- 
sively piled. In the multiplex model, however, between 
two successive receipts of grains, a node can topple due 
to an instability in the other layer. Accordingly, the like- 
lihood that a node becomes unstable as it receives a grain 
becomes smaller than l/Zc, and we numerically find that 
it scales roughly as 6>(A:i) ~ l/kf with a ^ 1.3 in our sim- 
ulation setting [Fig. 3(a), inset]. Therefore, the multiplex 
coupling induces truly nontrivial dynamics that cannot 
be viewed as a simple sum of two individual layers' dy- 
namics. 

The overall q^^'^ is shown in Fig. 3(a) along with the 
simplex counterpart, showing that in multiplex sandpile 
dynamics, branching into smaller number of branches 
(toppling a smaller number of grains) occurs more fre- 
quently than in simplex sandpiles. Meanwhile, the mean 
branching number 5 is 1 in both cases; numerically, 
we measured B^®'^ ^ 1.002 and B^^'^ ^ 1.005. This 
supports the robustness of self-organized criticality. Al- 
though ql^'^ takes a nontrivial form, its large-A: behavior 
is still dominated by the Poisson tail, which provides an 
explanation for the robustness of the avalanche size ex- 
ponent. 

To differentiate the cascade dynamics of the two mod- 
els in more detail, we measured the frequency of top- 
plings in nodes with ki and A:2, f{ki,k2) and plotted 
the normalized difference (j){ki^k2) — \f^®^{k\^k2) — 
/^®^(/ci, /c2)]//"^®^(^i, ^2) in Fig. 4. One can see that 
most nodes fail more frequently in the 102 model than 
in the 102 model whereas the nodes with links in only 
one layer (either /ci = or /c2 = 0) fail more in the 
simplex dynamics. This result suggests higher vulnera- 
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FIG. 5: (Color online) (a, b) Avalanche size distribution P{s) 
of the sandpiles on duplex networks with (a) two SF layers 
and (b) one SF and one ER layers. The degree exponent of 
the SF layers is 7 = 2.5. The slopes of both dotted lines are 
—5/3 and are drawn for the eye. (c, d) Branching probability 
Qk on (c) SF-SF and (d) SF-ER duplex networks. The slopes 
of both dotted lines are —2.5 and are drawn for the eye. 



bility of most nodes in multiplex sandpile cascades than 
in simplex scenarios. This picture is in line with recent 
studies indicating higher vulnerability of interdependent 
networks [lOj and multiplex cascade models [24] . 



B. On Scale- free Layers 

We have also performed numerical simulations on mul- 
tiplex networks with i) two scale-free (SF) layers and ii) 
one ER and one SF layer. For the SF layer, we used the 
static model [25 with the same mean degree z = 4 and 
the degree exponent 7 = 2.5. For the simplex dynam- 
ics, the degree distribution of aggregate networks in both 
cases has a power-law tail with exponent 7 = 2.5, which 
leads to the predicted avalanche size exponent r = 5/3 
according to Eq. (2). The multiplex dynamics is still 
critical (5^^^ ^ i ^102 ^ i.005), having the 

same power-law avalanche size distribution in both cases 
[Figs. 5(a,b)]. The detailed cascade dynamics are differ- 
ent; the branching probability qk is larger for very small 
and very large /c, but smaller for intermediate k in the 
multiplex dynamics than in simplex ones [Figs. 5(c,d)]. 
This result suggests an interesting picture that hubs be- 
comes more vulnerable in multiplex sandpiles. In an or- 
dinary (simplex) sandpile, hubs play the role of a shock- 
absorbers, by tolerating a large amount of stress (grains) 
without failing; in multiplex dynamics, hubs fail more 
often as there is an additional chance of failure due to 
its interdependency with a low degree node. Therefore, 
although the multiplex and the simplex dynamics do not 
differ in terms of scaling behaviors, the microscopic pic- 
ture of cascade dynamics at the individual node level 



- 10"' 

10-' 
10-^' 







; 


> 


AND-I A 
AND-II V 





(b)0 025 

0.02 
0.015 

0.01 
0.005 






m 








1(8)2 






102 □ " 






AND-I A 




g 


AND-II V 


1 







lO"" 10' 10^ 10^ 



10 15 

k 



FIG. 6: (Color online) (a) Avalanche size distribution P{s) of 
the AND-I (triangles) and AND-II (down triangles) models. 
The slope of the dotted line is —3/2 and is drawn for the eye. 
(b) Branching probability qk for the two models, compared 
with those of the 102 (circles) and the 102 (squares) models. 



does show differences that may have implications for risk 
assessment and aversion strategy for cascading failure in 
real-world phenomena. 



V. MODIFICATIONS OF THE MODEL 

In this section, we consider other variants of the BTW 
sandpile model on multiplex networks. The first variant 
is the one with a different rule iii) of the 102 model in 
Sec. III. In the new model, a node fails (topples) if it 
is unstable in both layers, that is, h£^x ^ for both 
i = 1,2. When a node fails, it topples C£^x = k^^x grains 
in each layer. We call this version the AND-I model. 
We also consider another rule such that a node topples 
all its grains. We call this variant the AND-II model. 
Therefore, for these variants we have instead of Eq. (5), 

hi^x ^ hi^x — C£,x (AND-I model), 

hi^x ^ (AND-II model). (8) 

For the AND-II model, when a node topples more than 
ki grains in the layer ^, it first distributes k£ grains to its 
neighbors and the remaining hi — k^ grains to randomly 
chosen neighbors. The scaling behavior is found to be 
still robust with respect to both variants of multiplex 
sandpile dynamics. They are critical (^^nd-i ^ 1.003 
and ^AND-II ^ 1.06) and have the same avalanche size 
exponent r = 3/2 on duplex ER networks [Fig. 6(a)]. 
Interestingly, the branching probability q^ for these two 
variants is found to be the same as that of the simplex 
1©2 model [Fig. 6(b)]. 

The final variant we consider is to change rule iv) of 
the 1(02 model such that a node fails if it is unstable in 
any layer, but the unstable node topples c^^x grains in 
both layers even when its current height is less than the 
threshold in layer 2. That is, we have, 

hi^x ^ rQax[/i£^a^ — Q,a?,0] (1(*)2 model) (9) 

for both and the height of all the neighbors of x in- 
creases by one. We symbolized this as the 1®2 model, 
to distinguish it from the original multiplex 102 model. 
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FIG. 7: Evolution of an infinite avalanche in the 1(*)2 model, 
(a) The total number of grains piled in the system at a spe- 
cific step of the infinite avalanche, (b) Number of unstable 
nodes at a specific step, (c) The accumulated number of nodes 
participating in the infinite avalanche up to the specific step, 
(d) The mean branching number B at the specific step of the 
infinite avalanche. 



In this version, grains are not conserved, but produced, 
during the avalanche. This setting is less physical, but 
can model situations in which the failure due to inter- 
dependency induces the full stress {c2^x) rather than the 
current stress level. We found that in this model, the sys- 
tem quickly enters into an infinite avalanche, in which the 
number of grains increases abruptly and eventually the 
avalanche covers the whole system, after which a finite 
fraction of nodes topple ad infinitum [Figs. 7(a-c)]. The 
branching analysis indeed shows that the mean branch- 
ing number becomes larger than one, increasing up to 
B ^ 4 [Fig. 7(d)]. Therefore, the branching process be- 
comes highly supercritical and can continue forever. The 
excess stress (grains) generated by interdependency re- 
generates further excess stress, and the system falls into 
the trap of an infinite sequence of cascading failure. The 
system is no longer critical. A similar absence of scaling 
had been observed for the conventional BTW model on 
Euclidean lattices without dissipation at boundaries [26] . 



VI. SUMMARY AND CONCLUSION 

We have studied a number of variants of the Bak-Tang- 
Wiesenfeld sandpile model on multiplex networks. To ac- 
count for the multiplexity of real- world complex systems, 
we introduced the multiplex sandpile model in which the 
network layers are interdependent so that a node fails if 
it becomes unstable in any layer (called the l(g)2 model). 
We found that this multiplex dynamics did not alter the 
scaling behavior of the mean-field-type sandpile dynam- 
ics, yet the detailed pattern of the microscopic cascade 
dynamics at the individual node level was affected by the 
multiplexity. Compared to the simplex counterpart, in 
multiplex dynamics, higher-degree nodes such as hubs in 
scale-free layers are more likely to fail. This finding res- 
onates with the recent finding of higher vulnerability of 
scale-free interdependent networks than Poisson-random 
ones under mutual percolation [10 . The robustness of 
scaling behaviors is further tested by considering other 
variants of multiplex dynamics. We could alter the scal- 
ing behavior only by breaking the conservation law, in- 
troducing excess stress due to interdependency, which 
drives the system supercritical. In this sense, the scal- 
ing properties of the BTW sandpile model is strongly 
robust on networks. In order to alter its scaling proper- 
ties, one needs to alter the asymptotic form of qk^ which 
is strongly constrained by the underlying substrate net- 
work structure. Therefore, random multiplex couplings 
between finite number of layers would not likely induce 
a nontrivial scaling behavior of BTW-type sandpile dy- 
namics despite differences in the microscopic dynamics. 
We note, however, that for other classes of cascade dy- 
namics, the multiplexity can play a critical role [ 22^ . 
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